arXiv: 1508.02584v2 [astro-ph.GA] 21 Aug 2015 


Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 24 August 2015 (MN [AT E X style file v2.2) 


Haloes light and dark: dynamical models of the stellar halo and 
constraints on the mass of the Galaxy 

A.A. Williams 1 *, N.W. Evans 1 

1 Institute ofAstronomy, University of Cambridge, Madingley Road, Cambridge, CB3 OH A, UK 


24 August 2015 


ABSTRACT 

We develop a flexible set of action-based distribution functions (DFs) for stellar haloes. The 
DFs have five free parameters, controlling the inner and outer density slope, break radius, 
flattening and anisotropy respectively. The DFs generate flattened stellar haloes with a rapidly 
varying logarithmic slope in density, as well as a spherically aligned velocity ellipsoid with a 
long axis that points towards the Galactic centre - all attributes possessed by the stellar halo 
of the Milky Way. We use our action-based distribution function to model the blue horizontal 
branch stars extracted from the Sloan Digital Sky Survey as stellar halo tracers in a spherical 
Galactic potential. As the selection function is hard to model, we fix the density law from 
earlier studies and solve for the anisotropy and gravitational potential parameters. Our best fit 
model has a velocity anisotropy that becomes more radially anisotropic on moving outwards. 
It changes from [3 ~ 0.4 at Galactocentric radius of 15 kpc to ~ 0.7 at 60 kpc. This is a 
gentler increase than is typically found in simulations of stellar haloes built from the mutiple 
accretion of smaller systems. We find the potential corresponds to an almost flat rotation 
curve with amplitude of ~ 200 kms 1 at these distances. This implies an enclosed mass of 
~ 4.5 x 10 n M 0 within a spherical shell of radius 50 kpc. 

Key words: Galaxy: halo - Galaxy: kinematics and dynamics - galaxies: kinematics and 
dynamics 


1 INTRODUCTION 


It has long been known that the distribution functions (DFs) of col¬ 
lisionless stellar systems depend on the integrals of motion. This 
result is often called |jeans| ( | 19 19} Theorem, though it was known 
earlier to Poincare. A number of choices of integrals of motion are 
possible. These include the classical integrals such as energy or 
angular momentum (e.g., [Eddington 1915 ; [Hunter & Qian|| 1993) 
| Evans & An|2006| l, the turning points or apocentric and pericentric 
distances of orbits l Hunter & de Zeeuw|1992[ >, or the actions (e.g., 
|Binney & Spergel|1984[ |B inne y|1987) . Although there are advan¬ 
tages and disadvantages to any of these choices, the actions remain 
very appealing because they are adiabatically invariant and they 
form a canonical set of coordinates when combined with the con¬ 
jugate angles (see e.g., |Goldstein|1980| . 

Historically, actions have not been much used in stellar dy¬ 
namics because of the difficulty in calculating them. For classi¬ 
cal problems like the Keplerian potential or the harmonic oscilla¬ 
tor, they can be readily evaluated by residue calculus (Bom|1927) . 
The isochrone is the most general spherical potential for which the 
actions are available as elementary functions (e.g., Henon 
[Evans et al.| 1990[>. For the separable or Stackel potentials, they are 
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available as a quadrature and provide a beautiful demarcation of 
phase space in terms of the orbital types ( |de Zeeuw|1985| ). How¬ 
ever, the last few years have seen a major change, in that a number 
of methods have been proposed to enable the rapid evaluation of ac¬ 
tions in a variety of potentials. This includes the adiabatic approxi- 
matio n (|Birmey|20l0 1 , local fitting (|San ders|2012]l and the Stackel 
fudge l |Binney|2012a |Sanders & Binney|2015a i. Together with the 
release of software packages like galpy i Bovy|2015] >, this work has 
transformed the calculation of actions into a routine matter for most 
axisymmetric and triaxial potentials. 


The form of the DF in action space is now an interesting prob¬ 
lem to study. The complex nature of galaxies suggests that each 
component - bulge, disk, stellar halo and dark halo — needs its own 
distribution function and is described by a characteristic shape in 
action space (Binne yj~1987[ . The functional form of the DF for thin 
and thick stellar disks is suggested by warming up models based on 
pure circular or epicyclic orbits and has been studied extensively 
in recent years, motivated by applications to the Milky Way (see 
e.g., |Binney|2010[rBovy & Rix|20 1 3[ . [Binney| i f20 1 2b[ showed that 
action-based DFs fit to the Geneva-Copenhagen Survey of nearby 
stars ([Nordstrom et al.|2004|l accurately predict the kinematics of 
stars in the deeper Radial velocity experiment (RAVE, Steinmetz] 
|et al.|2 006). Subsequently, Sanders & Binney (2015b) introduced 
extended DFs for the Galactic disk by introducing an analytic de- 
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dendence on metallicity. However, hot components such as bulges, 
dark matter haloes and stellar haloes have received much less atten¬ 
tion. |Posti et ah| (2015) and |Williams & Evans| 12015) considered 
self-consistent models with double power-law density profiles. In 
principle, such DFs can be used for stellar and dark haloes, as well 
as elliptical and dwarf spheroidal galaxies. However, these papers 
did not provide any fits to large datasets. We aim to remedy this de¬ 
ficiency here, by using such models to represent the Milky Way’s 
stellar halo as a tracer population in the Galactic potential. 

Section 2 introduces our family of DFs which depend on the 
actions. These are motivated by the density law of the stellar halo, 
which is of a flattened double power-law form to a reasonable ap¬ 
proximation. Section 3 provides a brief description of the data on 
blue horizontal branch stars extracted by |Xue et al.| ( |201 f) from the 
Sloan Digital Sky Survey. As the selection function of the sample 
is unknowable, Section 4 outlines our fitting methodology. Section 
5 presents our results, together with a comparison with other recent 
work on the Milky Way stellar halo. We sum up in Section 6, and 
suggest further extensions and applications of our models. 


2 THE DYNAMICAL MODEL 

In this section, we describe the construction of our dynamical 
model. First, we summarise the observational constraints already 
placed on the properties of the distribution function. We then 
demonstrate how these constraints can be used to construct a suit¬ 
able model. 


2.1 Observational motivation 


The spatial properties of the stellar halo are reasonably well- 
studied. |Deason et al.| (201lb) found that a sample of ~ 20000 
Blue Horizontal Branch (BHB) and Blue Straggler (BS) stars taken 
from the Sloan Digital Sky Survey (SDSS) is well modelled by a 
flattened broken power law 



p(r) oc 



r q < r b , 

r q > r b . 


0 ) 


where r 2 = R 2 + ( zlq ) 2 . Their best-fit parameters are r b = 27kpc, 
Ofi„ = 2.3, ttout = 4.6 and q = 0.59. Other studies have provided 
similar, though not identical, results (see e.g., |Pila-Dfez et al.|2015) 
and so we can be confident that the stellar-halo density is well- 
represented by oblate models with a sharp change in gradient in the 
density profile at ~ 30kpc. This rapidly changing gradient in the 
density is equally well described by the Einasto density law 


p(r q ) * exp |(r,/r eff ) 1/ " - l]j, (2) 

where d n = 3n - 1/3 + 0.0079/n. |Deason et al.|j201 lb) also fit such 
a profile in their analysis, with best fit parameters are ( n , q , r e g) = 
(1.7,0.58,20 kpc). 

The structure of halo stars in velocity space is not as well con¬ 
strained, which is primarily explained by the comparative lack of 
information: we possess all three spatial coordinates, but usually 
only one component of the velocity (the line-of-sight velocity), of 
each halo star. Nonetheless, the fact that the sun resides at a dis¬ 
tance from the Galactic Centre still allows us to place constraints 
on the velocity ellipsoid of the stellar halo.|Bond et al.|(|2010|> used 


a large sample of SDSS halo stars to investigate the velocity struc¬ 
ture within r ac ~ 10 kpc, and found that the velocity ellipsoid is 
essentially spatially invariant within their distance limit. Their in¬ 
ferred ellipsoid aligns with spherical coordinates, with axes given 
by cr r = 141 kms -1 , cr^ = 85 kms 1 and cr s = 75 kms -1 . The spher¬ 
ical anisotropy parameter, given by 


m = i - 


2 cr? 


(3) 


has a value of ~ 0.65. Other studies of halo kinematics find similar 
results, with p usually between 0.4 and 0.6 (e.g.,|Smith et al.|2009[l. 


2.2 Distribution function and gravitational potential 


Given the observational constraints, we wish to construct a DF that 
generates a flattened density profile with a rapidly changing log¬ 
arithmic slope, and with a velocity ellipsoid that is is spherically 
aligned and radially biased. Our model consists of a DF for the 
BHB population and a parameterisation of the Galactic gravita¬ 
tional potential. We choose to model the potential in a very simple 
way, using a spherical power law 


O(r) = -4 


S \ 10 kpc 


(4) 


where 0 < S < 1. Although the true Galactic potential is signifi¬ 
cantly more complex, this simple approximation is motivated by re¬ 
cent studies of the Galactic potential far from the disk. For example, 
in a recent study of the GD1 stream, |Bowden et al.[ (2015) found 
the potential to be well described by a scale-free and modestly flat¬ 
tened {q = 0.91 ± 0.04) model at distances r GC ~ 15 kpc. Using 
the same stream generation method, [Gibbons et al.|(20l4) showed 
that the mass profiles of flattened model galaxies are nonetheless 
accurately recovered using a model potential that is spherical. Fur¬ 
thermore, the spherical alignment of the velocity ellipsoid is also 
suggestive that the gravitational potential is itself close to spheri¬ 
cally symmetric at large radii (Evans et al. 2015). That said, our 
over-simplification of the Galactic potential is only permissible in 
a first attempt to fit data with action-based DFs, and should be su¬ 
perseded in later work. In particular, evidence on the shape of the 
dark matter halo of the Galaxy is confused, and trivial as well as 
axisymmetric shapes remain possible (e.gJLaw et al.[2009||Law &| 

|Majewski|20~r0l|Vera-Ciro & Helmi|20l3l(peg & Widrow|20 1 3) . 

In a spherical potential, the actions are the azimuthal action 
J#, which is equal to the angular momentum about the z-axis; the 
latitudinal action J g = L - \J^\, where L is the specific angular 
momentum, and the radial action 

J r=^f V2 [E - ©(r)] - ntf-Ar. (5) 

Since our choice of gravitational potential is a power-law, the result 
found by Williams et al. (2014) applies, which gives the Hamilto¬ 
nian of a generic power law as a function of the actions 


H(J) = C (L + DJ r y = C£\ (6) 


where C, D and £ are functions of the power law index 8 (for full 
expressions see Will iams et al.||2014) . |Williams & Evans| (2015) 
then showed that a model with two power-law regimes is generated 
by a DF of the form 


/(/)« 


£- 


{£ 2 + J b 2 ) 


2 \ 0 /-/ 0/2 ’ 


(7) 


where J b is the ‘break action’, and controls the location of the 
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break in the density profile of the model. We expect our DF to 
broadly resemble this ansatz, although here we are modelling a 
tracer population in an independent potential rather than building a 
self-consistent model. The DF of Equation Q generates a spheri¬ 
cally symmetric, isotropic double power-law model in the potential 
of Equation 0. since X is very nearly a function of energy alone. 
The parameters fj and v are related to the logarithmic slopes of the 
stellar halo density profile oro and a„ by 

A = ({a 0 /S-3/2), 

fi = ((.aJS-3/2). (8) 

In order to give the model an anisotropic, spherically aligned ve¬ 
locity ellipsoid, we follow Williams & Evans (2015), who demon¬ 
strated that modifying the argument X to the DF such that 

X —> L + f DJ (9) 
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Figure 1. Distributions of stars in the BHB dataset. Top: Histogram of 
Galactocentric distance. Middle: Distribution of Galactocentric distance 
and line-of-sight velocity. Note there is a noticeable overdensity at r ~ 
35 kpc due to Sagittarius stream members. Bottom: Positions on the sky of 
the BHB sample after our cuts have been made. The dashed horizontal lines 
are at |i>| = 30°, which define our cut in Galactic latitude in order to avoid 
disc stars. 


(/ > 0) endows the model with anisotropic kinematics. This arises 
because the factor f alters the relative importance of high and low 
eccentricity orbits. If 0 < / < 1, then the model becomes radially 
biased, and / > 1 gives tangentially biased models. 

Our final modifications to the DF will have the effect of break¬ 
ing the spherical symmetry of the model density profile, so that it 
becomes flattened. In action space, this means introducing an ex¬ 
plicit dependence on 7^, which will give the model a symmetry 
axis. |de Bruijne et al.| (: 1996) discovered a very elegant method for 
this purpose, though only for scale-free models. Remarkably, they 
found that introducing a multiplicative factor has precisely this ef¬ 
fect. This factor is given by 


h(e 2 jj/L 2 ) 


y dMf)* (e 2 J 2 A 

h *«& L1 ) ’ 


( 10 ) 


where e = -^1 — q 2 is the ellipticity, (...)*, is the Pochammer sym¬ 
bol and a is the power law slope in the density. This factor has 
the effect of reducing the number of stars on orbits where the ratio 
J<l>/(\Jip\ + Jo) is small, which reduces the amount of vertical mo¬ 
tion in the model relative to the circulation around the symmetry 
axis, thereby flattening the density profile. The expression at first 
sight appears unwieldy due to the sum to infinity, but numerical ex¬ 
periments demonstrate that the series converges very quickly, after 
~ 10 terms. As it stands, this expression is inappropriate for our 
purposes, since it is designed for scale-free densities. We make the 
replacement 


a —> a(J) - 


a 0 + aU\J\/Jb) 2 
1 + (\J\/Jb) 2 


(ID 


where \J\ = y 2 Jr. This approximately preserves the constant flat¬ 
tening of the model. A suitable DF is therefore 

N h(J)X- A 


fU) = 


where At is a normalisation factor such that 


( 12 ) 


(2?r) 3 J f(J)d 3 J = 1, (13) 

and h(J) is given by Equations © and We have now con¬ 
structed a model with all of the desired properties motivated by 
observations. All of our parameters are nicely physical, except for 
the break action. In practise, we should like to replace this param¬ 
eter with a break radius, r b . Unfortunately, a wholly analytical ex¬ 
pression for J b in terms of a break radius is unavailable. We can 
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Figure 2. Likelihood contours from our analysis. We show the 68% and 
95% confidence intervals of the marginalised distribution over pairs of 
model parameters, and mark the maximum likelihood gridpoint with a red 
cross. In the one-dimensional distributions, the maximum likelihood value 
for each parameter is marked with a dashed vertical line. One can see that 
the circular velocity at lOkpc is tightly constrained, and that models close 
to logarithmic potentials are favoured. The contours of marginalised dis¬ 
tributions involving 8 are one-sided due to the fact that we only consider 
declining rotation curves. Our grid in 8 was truncted at 8 = 10 -2 because of 
the highly singular limit that occurs at 8 = 0. 



approximate J h with the expression 

J b = ( vo \ r £ _ -ul 1/f ( r o 

lOOOkms -1 kpc llOOkms^ 1 / yfg l 2 \ \ lOkpc 

(14) 

This expression is derived by considering the relative contributions 
of perfectly radial and circular orbits to the density profile of the 
model. For a given r 0 and Vo, the model density profile will break 
at an elliptical radius r\, - ro if Jb from the above expression is 
used as the break action. We then choose to sample the parameter 
r 0 . The model parameters are then 

V = (ffo, a m , ro, q, /, v 0 , S ), (15) 

each of which has a simple physical interpretation. 


3 DATA 

In this study, we fit our model to the dataset provided by |Xue et al.| 
< ]20 1 1} , who compiled a catalogue of BHB stars from the SDSS 
DR8 catalogue ( |Aihara et al,|20 1 1 The dataset contains ~ 4000 
spectroscopically confirmed BF1B candidates, arising from various 
spectroscopic surveys within SDSS. BHB stars are a very useful 
class of object to use when investigating the stellar halo, as they 
are intrinsically bright, with M g ~ 0.5, and Deason et al.| l |2011b) 
showed that very accurate estimates of their absolute magnitudes 
can be derived (A M g — 0.15) via the polynomial relation 

M g = 0.434 - 0.169(g — r) + 2.319(g - r) 2 

+ 20.449(g - r) 3 + 94.517(g - r) 4 , (16) 



dlos (kms ') 


Figure 4. Comparisons of our maximum-likelihood model to the data. Top 
panel: the line-of-sight velocity dispersions of the data and mock catalogue 
as a function of Galactocentric radius. The two profiles are in excellent 
agreement, other than an upturn in the profile of the data at large radii. Mid¬ 
dle: residuals between the model and the data in the dsitribution of Galacto¬ 
centric radius and line-of-sight velocity. The residuals are essentially noise 

other than a feature between 30 and 50kpc at-100 kms -1 , which is due 

to stars belonging to the leading arm of the Sagittarius stream. Bottom: 
comparison of the distributions in line-of-sight velocity between the model 
and the data. The two distributions are in good agreement, except the data 
displays a slight skew towards negative line-of-sight velocities, due to the 
overdensity of Sagittarius stars at vlos -100 kms -1 . 
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Figure 3. Spatial properties of the model. Left: the density profile of the model in the x—z plane. Right: the ID density profile of the model with elliptical 
radius. The black dashed line denotes the value of the parameter ro, and the grey lines surround the region for which we have data. Overplotted in the range 
where the data lies is the best-fit Einasto profile from |Deason et al. i foTThj . which we use to constrain the density profile of our model. 


which is valid in the colour range -0.25 < g - r < 0.0. From the 
absolute magnitude, we can estimate distances to the stars in our 
sample to high precision. Before carrying out our analysis, we make 
two cuts on the dataset. First, we limit the colours to the range for 
which Equation is valid, to ensure that the inferred distances 
are reliable. We then only include stars for which \b\ > 30° in order 
to reduce contamination from the disk. This leaves us with 3534 
stars at Galactocentric radii 10 kpc < r < 50 kpc. each of which has 
4 of 6 phase-space coordinates very accurately constrained. Figure 
[TJdepicts the distributions in Galactocentric distance, line-of-sight 
velocity and on-sky positions of the data. 


4 FITTING METHODOLOGY 

Ideally, we would like to use this dataset to constrain all the pa¬ 
rameters of our model. This would mean inferring properties of 
the entire phase-space distribution of the stellar halo, including the 
density profile and kinematics. However, in order to carry out such 
a study, we would require detailed knowledge of the selection func¬ 
tion of the sample of stars that is being used. The observed distri¬ 
bution of stars in phase-space and chemistry, Z, is the product of 
the selection function with the true density, i.e. 

/ obs (x, v,Z) = S (m, l, b, Z)x f mc (x, v,Z). (17) 


then proceed to use the velocities of the stars to constrain the kine¬ 
matics of the stellar halo, as well as the Galactic mass profile. This 
is precisely what was done by |Deason et al.|(2012| hereafter D12), 
though they used a smaller subsample of ~ 1900 stars from this 
dataset. Our DF lets us model the entire dataset, and thus place 
tighter constraints on halo kinematics and the mass of the Galaxy, 
though of course the properties of the stellar halo such as the flat¬ 
tening (q = 0.59) are fixed at outset. 

Given the above discussion, we use the data to constrain only 3 
of the 7 model parameters: /, which controls the kinematics, Vo, the 
circular speed at 10 kpc, and 6, the slope of the gravitational poten¬ 
tial. The remaining parameters are fixed for each set of (/, t’o , S) 
by optimizing the fit of our model’s density profile to the Einasto 
model of Deason et al. (201 lb) described in Section [2~T| This stage 
of the procedure is equivalent to enforcing a very strong prior on 
the remaining 4 parameters of the model. The likelihood for an in¬ 
dividual datum is then given by 


P(VLOS,il*i> /> v 0 , tf) = 


/ 


//</> 

//</> 


d(j, e) 
d(X) 
d(J, 0) 
d(X) 


dv ; di’* 


dv, dv*dv LO s, 


(18) 


where X are the usual Galactic coordinates. The above expression 
simplifies, owing to the fact that the Jacobian factors depend only 
on position coordinates, giving 


Note that the selection function, S , is presumed not to depend on 
the velocities of the stars, but does depend on their magnitudes (m), 
on-sky positions and chemistry. Without an effective model of S , 
we cannot hope to constrain the spatial and chemical distributions 
of the stars. Unfortunately, the selection function for our data is 
essentially unknown. Since the stars were observed on a variety 
of different plates from different spectroscopic surveys, and they 
have a relatively low on-sky density, we concluded that a reverse- 
engineered selection function of the sort used by |Bovy et al. | i |20 Y2\ 
cannot be produced. This means inference on the density profile of 
the stellar halo is not possible with the data used here. 

All hope is not lost, however. Although some of the parameters 
of our model cannot be constrained by this data, we can still pro¬ 
ceed. We can use the independent analysis of |Deason et al.|(20 1 1 bfr 
to place strong priors on the spatial distribution of the stars, and 


P(vlos,i I Xi , /, v 0 , 5) = —— f f(J) dv, dv*, (19) 
p(Xi) J 

which is the normalised line-of-sight velocity distribution at the 
position of the star in question. The log-likelihood for the entire 
dataset is then simply 

Nbhb 

log/ = Xj lo SPj- (20) 

7=1 

In order to compute this likelihood with minimal noise, high- 
precision numerics are required for a number of different tasks. 
First, we need to evaluate the non-classical integral J r in order to 
compute the integrand /(/). In order to do this, we compute J, nu¬ 
merically using adaptive Gaussian quadrature implemented in the 
Gnu scientific library over a grid of 200 points in angular momen- 
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turn at some fiducial energy E 0 . Since the potential is scale-free, the 
radial action scales self-similarly with energy, so we do not need to 
consider other energies: J r can be recovered by simple multiplica¬ 
tion by a scale factor. In-between grid points, we recover J r using 
cubic spline interpolation. 

Once the radial action is available, we must perform the two 
(marginalising over the proper motions) and three-dimensional (in¬ 
tegrating over velocities to find the density) integrals required to 
compute the normalised line-of-sight velocity distributions. This 
is done using the Cuhre algorithm implemented in the Cuba pack¬ 
age QHahn|2005] >. This particular algorithm is deterministic, which 
guarantees that our likelihoods are replicable. 

Evaluating the triple integral needed to calculate the density of 
the model at a given position is significantly more time-consuming 
than the double integral when marginalising over the proper mo¬ 
tions. However, the model has a smooth two-dimensional density 
profile, and we can use this to our advantage. Instead of explicitly 
evaluating the density at the positions of all 3534 stars, we instead 
construct a grid and use bilinear interpolation to recover the density. 
An astute choice of variables proves to be 

u = ilog(fl 2 + (zA?) 2 ), (21) 

(qR\ 

t = arctan — . 

Ui/ 

An economically sized grid can be used for interpolation to find the 
density. We use a grid of 20 points in u and 10 points in f, and find 
that the density at the position of any star is recovered to a precision 
of 0.4 percent at worst. 

All of the above calculations are performed in C, which is 
then wrapped for use in python. In order to locate the maximum- 
likelihood solution, we perform a grid-search. Monte-Carlo meth¬ 
ods are unnecessary in this instance, due to the low dimensionality 
of the parameter space. 


5 RESULTS & DISCUSSION 

Here, we present the results of our analysis. First, we consider how 
well the model fits the data. Then we discuss the constraints placed 
on the Milky Way mass distribution and circular velocity profile, 
comparing our results with those from other analyses, finally, we 
discuss the kinematics of our best-fit model. 

5.1 Comparison to the data 

Figurepldepicts the likelihood distributions from our analysis, and 
Table [nsummarises the best-fit parameters and their errors. Note 
that four of the parameters are fixed by our prior on the density pro¬ 
file, and Figure [3] depicts the one and two-dimensional properties 
of the density of this model. 

In order to compare the model to the data, we generated a 
mock catalogue by drawing a line-of-sight velocity at the position 
in the Galaxy of each star from the maximum likelihood model. 
The reason we do not also sample position is because, as previ¬ 
ously stated, the selection function is unknown and depends on the 
location of the star. We have assumed that the selection function 
does not depend on the kinematics of a star, and we are therefore 
safe to freely sample the line-of-sight velocity distributions of the 
model at positions where stars have been observed. We draw the 
mock catalogue using a simple rejection sampling technique. 

Figure [4] depicts several comparisons between the data and 


our mock catalogue. The upper panel, depicting the line-of-sight 
velocity dispersions as a function of galactocentric radius, shows 
excellent agreement between our model and the data. There is an 
upturn in the profile of the data at large radii, which is due to the 
presence of stars belonging to the Sagittarius stream. This sub¬ 
structure should then be apparent in the two-dimensional distri¬ 
bution in galactocentric radius and line-of-sight velocity. Indeed, 
as the middle panel of Figure [^demonstrates, there is a signal at 

the same distances at vlos -lOOkms -1 . The final panel depicts 

the overall distribution in line-of-sight velocity. The profiles are in 
good agreement, except the data shows a slight skew towards neg¬ 
ative line-of-sight velocities, which is again due to the presence of 
Sagittarius members. 

We confirmed that the systematic residuals between the data 
and the mock catalogue are a consequence of the Sagittarius stream 
by making a further cut to the dataset. We removed stars with \B\ < 
10°, where B is one of the two Sagittarius stream coordinates (A, B) 
defined in Belokurov et al. (2014). First, we inspected the line- 
of-sight velocity distribution of the stars beyond r = 30 kpc, and 
found that the skew towards negative radial velocities was removed 
by this cut. We then drew another mock catalogue from the model 
and compared it to the reduced dataset. The signal seen at Vlos ~ 
-100 kms -1 is no longer present, confirming that the contamination 
in the sample is largely from the Sagittarius stream. 

Overall, the model seems to be in good agreement with the 
data, and we are confident that the Sagittarius stream members do 
not bias our fits, since the majority of our constraining power comes 
from stars at galactocentric radii < 30 kpc, for which we do not 
possess significant contamination. In fact, this highlights an inter¬ 
esting application of smooth models of the stellar halo: one can fit 
a smooth model to the data, and search for substructure using the 
residuals between the data and the model. Having established that 
our model is as good representation of the data, we can now go on 
to discuss the implications of our inference. 

5.2 The Cumulative Mass Profile of the Milky Way 

The parameters v 0 , the spherically averaged circular speed at 
10 kpc, and 8 , the slope of the potential, are tightly constrained by 
our analysis. v 0 is found to be ~ 200 kms -1 , and 8 is very close to 
zero, implying the rotation curve is almost flat. Taking these results 
together, we find that the predicted enclosed mass within a spheri¬ 
cal shell of radius 50kpc is 4.48+jj x 1O U M 0 . 

Figure[5]depicts the cumulative mass distribution and circular 
velocity curves of the Milky Way predicted by our model, along 
with the 68% confidence intervals. We have also plotted the results 
from other studies on this matter, which we shall now discuss in 
detail. There are multiple predictions of the Milky Way enclosed 
mass at ~ 50kpc. D12 and |Xue et al.| (2008] hereafter X08) both 
also use samples of BHB stars from the SDSS: D12 used a sub¬ 
sample of ~ 1900 stars from the same dataset considered here, and 
X08 used a sample from an earlier data release. Our result agrees 
extremely well with that of D12, who also used distribution func¬ 
tion based modelling to place constraints on the mass profile. As 
one might expect, our error-bar is somewhat smaller, owing to the 
fact that we use roughly double the amount of data in our study. 
Interestingly, the logarithmic slope of the potential inferred by D12 
is larger than our result: 8 ~ 0.4. Their value is almost half-way 
between the Keplerian and logarithmic cases, so the rotation curve 
has a steeper decline with radius. This is an interesting result, and 
has several possible explanations. It is clear that the gravitational 
potential of the Galaxy cannot possibly be as simple as is implied 
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Table 1. Parameters of our maximum likelihood model and their dispersions. We do not quote a dispersion for the flattening, q, because it is kept fixed during 
the fitting procedure. 




Figure 5. Maximum likelihood profiles (red curves) and 68% confidence intervals (grey shade) for the mass profile (left) and circular velocity curve (right) of 
the Milky Way. The error bars come from sampling the 1 cr confidence interval on the joint posterior of the two potential parameters. Dashed lines surround 
the region in Galactocentric radius for which we possess data. Also plotted are the quoted error bars from several other studies of the Milky Way cumulative 
mass distribution: D12 (Deason et al.|2012) , G14 (Gibbons et al.|2014) , W99 fWilkinson & Evans|1999) , X08 (Xue et al.|2008) , G10 (Gnedin et al.|2010) and 
W10 i Watkins et al. 2010J. The markers corresponding to D12, G14 and W99 are all offset from 50kpc so that they are distinguishable. 


by Equation Q. and so biases could certainly be introduced into 
an analysis that assumes such a parameterisation. In particular, if 
the slope of the Galactic rotation curve changes appreciably over 
the range of distances between 10 and 50kpc, the lack of flexibility 
in the model means this cannot be accounted for. D12 only consid¬ 
ered stars beyond an elliptical radius r q = 27 kpc, the more distant 
stars in this dataset, so it is not unreasonable to assume that their 
analysis was more sensitive to the slope of the potential at larger 
distances. In our case, the left-hand panel of Figure ^informs us 
that the majority of the stars reside at radii 10 kpc < r < 25 kpc. 
Therefore, it may possibly be the case that the rotation curve is ap¬ 
proximately flat between 10 and 25 kpc, but begins to decline more 
sharply thereafter. 

X08 took a rather different approach in their modelling of the 
BHB population. Rather than fitting analytical distribution func¬ 
tions, they instead compared their data to mock samples drawn 
from cosmological simulations. Although they also predict a very 
flat rotation curve, its amplitude is somewhat lower. This discrep¬ 
ancy could be accounted for by the above reasoning, where the cir¬ 
cular velocity begins to decline more rapidly with radius. 

Wilkinson & Evans (1999, hereafter W99) estimated the mass 
enclosed at 50 kpc by modelling the distribution of the satellite 
galaxies and globular clusters of the Milky Way. They used con¬ 
stant anisotropy tracer distribution functions embedded in spher¬ 
ically symmetric models with rotation curves that are flat in the 
inner parts, then decline in a Keplerian fashion at large radii. These 
models have a circular velocity of form: 

, vl 

v? = . ° (22) 

yl + r 2 /a 2 

Their result, M(r < 50kpc) = 5.4^ x 10 u M o , is in excellent 


agreement with ours (note the asymmetric error-bars). It is also 
interesting to note that the scale-radius of their model is found to 
be very large (140 kpc < a < 260 kpc), implying a very flat rotation 
curve at the distances of the stars in our dataset. 

|Gibbons et al.| ( f2014| > used observations of the Sagittarius 
stream, in particular the precession of the apocentric position of 
the stream, to infer the mass profile between radii 50 kpc < r < 
100 kpc. As the title of their paper states, they find a particu¬ 
larly ‘skinny’ Milky Way, with an enclosed mass at 50 kpc of just 
2.9±0.4x 10 n M o . Given the small error bars on their measurement, 
their result is in tension with ours. Evidently, the two methods pos¬ 
sess different systematic biases. 

We now compare results from two other studies, [Gnedin et aT| 
( |2010| hereafter G10) and |Watkins et al.| ( |2010[ hereafter W10), 
who both inferred the mass enclosed at larger distances. Although 
it is still interesting to discuss these results, our result must be ex¬ 
trapolated in order for comparisons to be made. G10 used a large 
sample of radial velocity measurements of halo stars to carry out 
a Jeans analysis. They assume a spherically symmetric tracer den¬ 
sity with logarithmic slope between -3.5 and -4.5 with a constant 
anisotropy, and find M(r < 80kpc) = 6.9^® x 10 u M o . In spite of 
their less sophisticated model for the tracer density, their result is 
in good agreement with ours. Finally, Watk ins et al.| ( [20T0l l took a 
sample of 26 satellite galaxies of the Milky Way and applied their 
virial mass estimators to the data. The estimators assume a power- 
law tracer density, as well as a power-law gravitational potential. 
Their measure is very sensitive to the assumed anisotropy of the 
tracer population, and so their estimate for the enclosed mass has 
a large error bar. Nonetheless, we are consistent with their conclu¬ 
sions at the 68% confidence level. 
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Figure 6. Velocity ellipsoids of the best-fit model across a variety of po¬ 
sitions. The ellipsoids are everywhere aligned with spherical coordinates, 
with long axes directed towards the Galactic centre. The size of the ellip¬ 
soids do not vary a great deal with radius, suggesting a relatively flat ve¬ 
locity dispersion profile. Their shape becomes more elongated with radius, 
indicating that the velocity anisotropy of the BHB population is growing 
with distance from the Galactic centre. 


5.3 Kinematics of the Stellar Halo 

We now study the kinematics of our best-fit model. Figure [^de¬ 
picts the velocity ellipsoids of our model over a range of positions 
in the Galaxy. Our model naturally provides a velocity ellipsoid 
that is everywhere aligned with spherical coordinates, in agreement 
with the observations (e.g. |Smith et~aT|2009| |Bond et ak||2010| >. 
The ellipsoid always has a long-axis pointing towards the Galactic 
centre, meaning that the orbital distribution is everywhere radially 
biased. The ellipsoid becomes more elongated with distance from 
the Galactic centre, implying that the stellar halo becomes more ra¬ 
dially biased in the outskirts. Each ellipse is coloured by the value 
of the anisotropy parameter 

a 2 . + a\ 

m = i - 0-0, (23) 

2o-2 

which is a measure of the local velocity anisotropy. If /? > 0 (< 0), 
then the model is radially (tangentially) biased, and /3 = 0 implies 
an isotropic model. Figure[7]depicts maps of the three spherical ve¬ 
locity dispersions with position. We can see that the radial and az¬ 
imuthal dispersions are oblate, whereas the latitudinal dispersions 
are prolate in their distribution. The latitudinal dispersions also ex¬ 
hibit a ‘pinching’ in their profile in the equatorial plane. This fea¬ 
ture appears generic in radially biased models that are oblate or 
triaxial (San ders & Ev ansj2 0l5|l. 

An obvious comparison to make is between the prediction of 
D12, whose DF enforced a globally constant anisotropy, and that 
of our model. Figure [8] depicts the anisotropy parameter and its 
uncertainty as a function of distance along the direction R = z in our 
model, as well as the value picked out by D12 and their confidence 
intervals. The two profiles are in agreement within the error bars 
at all positions for which we have data. Note too that D12’s choice 
of DFs was questioned by |Fermani & Schdnrich||2013j >, but our 
analysis seems to vindicate the work. 

Our model predicts an increasing anisotropy profile with ra¬ 
dius, but given the sparsity of the current data, it is difficult to as¬ 
sess how realistic this is, so we turn to numerical simulations to 
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Figure 7. Maps of the velocity dispersions cr r (top), <r<s (middle) and erg 
(bottom). The radial and azimuthal velocity dispersion contours are oblate, 
whereas the latitudinal velocity dispersions are prolate. The radial velocity 
dispersions decrease the least with Galactocentric distance, and the latitudi¬ 
nal dispersions the most — as we expect for a flattened density. 


make a comparison. Bullock & Johnston (20051 simulated 11 dif¬ 
ferent stellar haloes in the ACDM context, each with the same dark 
halo mass and stellar disk at redshift zero, but with distinct accre¬ 
tion histories. The stellar haloes are built up over cosmic time by 
the accumulation of many different subhaloes which have had stars 


‘painted’ onto them. 

Kafle et al. ( 2012 ) then used the Galaxia code 

(Sharma et al. 2011 

i to construct synthetic BHB populations for 


these simulations, and analysed their velocity anisotropy profiles. 
They find that the mean anisotropy profile of the 11 different BHB 
populations is well represented by a function 


m = 


Par 2 

rl + r V 


(24) 


where j3 0 = 0.765 and r a = 2.4kpc. A rising anisotropy profile is 
therefore consistent with their findings, although there are signif¬ 
icant differences between the simulation properties and those we 
infer here. We instead find that the anisotropy profile of our model 
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is well represented by the function 


J3(r)=j3o + 


(fli — Po)r 
r + r 0 


(25) 


withySo = 0.05, = 0.82 and r 0 = 18.2kpc. The anisotropy profile 

rises at a much lower rate with radius in our model than in the 
simulations, only reaching /? ~ 0.7 at r ~ 100 kpc. For reference, 
we also plot the anisotropy profile from the simulations in Figure[8] 
The differences are curious, since our model is capable of creating 
profiles similar to that of Equation (|24b, but they do not seem to 
be preferred. The Bullock & Johnston q2005 ( i simulations involve 
many small haloes accreting across cosmic time to build up the 
stellar halo, but this is by no means the only possibility for the 
formation history of the Galaxy’s halo. 

In the future, surveys such as Gaia will provide the three di¬ 
mensional velocity distributions of halo stars. In fact, such dat- 
acubes are already beginning to be available. |Bond et~ak] ( |2010| 
hereafter B10) analysed the motion of nearby stars from the SDSS 
using radial velocities coupled with proper motions derived from 
SDSS and Palomar Observatory Sky Survey (POSS) astrometry. 
Figure [9] depicts the two-dimensional velocity distributions of a 
mock catalogue of 2000 stars drawn from the best-fit model at the 
position R = 8 kpc, z = 3.5 kpc. We draw the mock sample using a 
rejection sampling method, with a sampling distribution that is the 
product of three broad Gaussians in velocity space. The distribu¬ 
tions are qualitatively very similar to those found by B10, although 
there are some differences (see their Figure 12). The halo distri¬ 
butions presented by B10 are broader in v z , which is somewhat 
surprising given that the halo is a flattened system. Furthermore, 
B10 infer a greater velocity anisotropy in this region,/! ~ 0.65, as 
compared with our model, which has /3 ~ 0.4. This is perhaps sug¬ 
gestive that the distribution of halo stars is significantly affected by 
the presence of the disk, which is not accounted for in our model. 
However, we have not convolved the velocities here with errors, or 
considered any contamination from disk stars, both of which will 
have an effect on the observed distributions. 


6 CONCLUSIONS 


We have presented a new model for the Milky Way stellar halo, 
and used it to fit a sample of Blue Horizontal Branch stars from the 
Sloan Digital Sky Survey i jXue et al.|201 l[ l. For the first time, we 
have used a distribution function formulated in terms of the action 
integrals in the fitting analysis. The density profile generated by 
this DF is flexible, which allows us to use the results on the SDSS 
photometry from Deason et al.| ( [201 lb| > as a prior on the spatial 
distribution of our BHB dataset. We argue that such a prior is nec¬ 
essary due to the unknown selection function of the data. We then 
fit three model parameters using the sample of 3534 stars, two cor¬ 
responding to a simple power-law gravitational potential and one 
that controls the kinematics of the model. 

Our results are consistent with a very flat rotation curve for the 
Milky Way galaxy, with a mass enclosed at 50kpc of 4.5x 1O U M 0 . 
This is in good agreement with other recent estimates in the lit¬ 
erature. The spherically aligned velocity ellipsoid of our model is 
everywhere radially elongated, with a radial bias that increases with 
Galactocentric distance as 


m = a > + 


(J3 1 -J3 0 )r 
r + r 0 


(26) 


with /3o = 0.05, = 0.82 and r 0 = 18.2 kpc. This behaviour qual¬ 

itatively agrees with results seen in numerical simulations, in that 



Figure 8. Comparison between the anisotropy predicted by our model 
(black line, grey shade in 68% confidence intervals), the best-fit value from 
D12 (dashed red line, red shade in 68% confidence intervals) and the Bul¬ 
lock & Johnston simulations (green line). Our more sophisticated model is 
in agreement with the analysis of D12 at the 68% confidence level, but in¬ 
consistent with the simulated haloes. The simulations are far more radially 
biased than our model, implying that there may be significant differences 
between the Milky Way stellar halo formation history and the typical pic¬ 
ture assumed in the ACDM paradigm. 


the anisotropy parameter is an increasing function of radius, but our 
model suggests a much gentler growth of /3 than in the simulations. 
Given that the model can produce far more radially biased distribu¬ 
tions than our maximum-likelihood model, it is interesting that the 
data do not seem to be consistent with the canonical simulations of 
stellar halo formation. This is suggestive that the formation history 
of the stellar halo may be somewhat different from that implied by 
the models of Bulloc k~& Johnston| ( |2005^ . 

Though our model is the most complex DF yet fitted to data 
on the stellar halo, it still has some obvious limitations. We have 
tightly constrained the density profile of the stars in our model so 
that it is in agreement with past analysis of the stellar halo, a step 
we unfortunately believe to be necessary owing to the unknown se¬ 
lection function of the data (c.f. Das & Binney, in prep.). Evidently, 
in applications to datasets provided by Milky Way surveys such as 
Gaia, we will have a better understanding of the selection function 
and be able to constrain the spatial and kinematic structure of the 
stars, as well as the Galactic potential. 

Another layer of sophistication can be added to our models 
in order to account for rotation in the stellar halo. Here, we have 
presumed that there is no net rotation of the halo, when in reality 
this score is by no means settled (e.g.,|Deason et al.|2011a||Fer-| 
|mani & Schonrich 20131. One way of introducing a mean azimuthal 
streaming velocity is by modifying the DF to become (e.g.,[Deason] 
|et al.|201 la||Binney|2014| > 

fUJ ) = [ytanh + (1 - y)] /(/), (27) 

where SJ is a small action to avoid discontinuities in the gradient 
of the DF and /(/) is the original ansatz from Equation ( | 1 2| ). The 
parameter y then takes values between -1/2 (maximal retrograde 
motion) and 1/2 (maximal prograde rotation). We can then add the 
extra parameter y to the analysis and search for rotation in the stel¬ 
lar halo. 

The model presented here is constructed in a very simple po¬ 
tential, meaning that it is possible to write down an ansatz for a DF 
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Figure 9. The local cylindrical velocity distributions from a mock catalogue 
at the location R = 8kpc, z = 3.5 kpc. We use a simple rejection sampling 
method to generate a catalogue from our model. Note the tilt of the veloc¬ 
ity distribution in the - v z projection, which results from the spherical 
alignment of the velocity ellipsoid. 


with predictable features. However, given the rapid recent devel¬ 
opments on action estimation in generic potentials (Binneyj2012a| 
|Bovy||20T4] |Sanders & Binney||2014| |Sanders & Binney||2015a| >, 

it is now of importance that we develop flexible DFs with well- 
understood properties in more realistic Galactic potentials. Current 
models of spheroidal components of galaxies ^Williams & Evans] 
|2015||Posti et al.|2 015 1 have been developed to produce certain fea¬ 
tures in spherical potentials, and even though they can be modified 
for use in axisymmetric potentials (e.g. Das & Binney, in prep.), 
their properties are less well understood in this case. Eventually, 
it will be the case that we simultaneously fit multiple-component 


Galaxy models (e.g. |Piffl et ah]|2015[ > to data from various sur¬ 
veys, but for this to be an effective procedure we must ensure that 
the DFs used for the various components are sufficiently sophisti¬ 
cated to model the high-quality data that is to appear over the next 
few years (Perryman et al.|200l) . For example, |Sanders & Evans 
1 201 5} in vestigate triaxial generalisations of the models in |Williams 
& Evans (20151, which should provide valuable information as to 
how to progress can be made in this area. 
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